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DIAGNOSIS  OF  DAMAGE  IN  SDF  STRUCTURES 


1.0  Introduction 

This  report  summarizes  the  results  of  an  investigation  in 
the  diagnosis  of  structural  damage. 

The  reliability  of  a  structure  is  the  probability  that  it 
will  survive  and  respond  to  inputs  in  a  satisfactory  manner  over 
a  preestablished  period  of  time.  If  the  time  through  which  a 
structure  is  intended  to  survive  is  chosen,  then  the  reliability 
of  a  structure  can  be  re-estimated  at  various  times;  and  the  re¬ 
liability  may  vary  as  a  function  of  time.  At  any  time,  the  reli¬ 
ability  of  a  structure  is  dependent  upon  the  class  of  inputs  to 
which  the  structure  will  be  subjected  and  the  condition  of  the 
structure.  If  damage  in  a  structure  accumulates  at  a  faster  rate 
than  anticipated  during  the  initial  reliability  estimation,  then 
its  reliability  will  diminish  with  time.  Repairs  can  be  per¬ 
formed  on  a  damaged  structure  to  increase  its  reliability. 

In  view  of  the  potential  variability  of  structural  reliabil¬ 
ity  and  the  dependence  of  reliability  on  the  damage  level  in  a 
structure,  it  is  desirable  to  know  how  much  damage  has  accumu¬ 
lated  in  a  structure  at  a  given  time.  When  damage  occurs  in  a 
structure,  it  can  appear  in  the  form  of  cracks  and  permanent 
structural  deformation. 

An  indicator  of  the  damage  which  accumulates  in  a  structure 
as  the  result  of  structural  motion  over  a  period  of  time  is  the 
energy  dissipated  by  the  structure  during  that  time.  Therefore, 
we  wish  to  estimate  this  quantity.  Particularly,  it  is  important 
to  estimate  the  energy  dissipated  in  a  structure  as  it  responds 
to  strong  motion  excitations,  for  this  type  of  response  causes 
the  greatest  damage.  The  energy  dissipated  by  a  structure  re¬ 
sponding  to  strong  motion  excitation  is  dissipated  due  to  the 
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hysteretic  behavior  of  the  structural  material.  The  equations 
governing  the  hysteretic  structural  response  of  a  lumped  mass 
system  are  second-order,  nonlinear,  ordinary  differential  equa¬ 
tions  with  hi  story -dependent  stiffness  terms. 

In  the  present  investigation,  it  is  assumed  that  the  hyster¬ 
etic  response  of  a  lumped  mass  system  is  approximately  governed 
by  a  higher-order,  linear,  ordinary  differential  equation.  This 
assumption  is  motivated  by  studies  summarized  in  the  literature. 
For  example.  Wen  [1]*  and  Baber  and  Wen  [2,3]  have  used  a  third- 
order,  linear,  ordinary  differential  equation  to  approximately 
govern  the  response  of  a  hysteretic  single-degree-of-f reedom 
(SDF)  system  to  white  noise  and  non-white  stationary  random  in¬ 
puts;  and  they  have  used  a  system  of  third-order,  linear,  differ¬ 
ential  equations  to  approximately  govern  the  response  of  hyster¬ 
etic  mul ti-degree-of-f reedom  (MDF)  structures.  Their  results 
show  that  the  higher-order  linear  system  provides  an  acceptable 
model  for  the  hysteretic  system.  Lutes  [4]  used  linear,  third- 
and  fourth-order  equations  to  approximately  govern  the  response 
of  an  SDF  hysteretic  system.  His  results  show  that  the  spectral 
density  of  the  hysteretic  structural  response  can  be  matched  by 
the  linear  systems;  in  addition,  the  buildup  of  structural  re¬ 
sponse  in  the  hysteretic  system  can  be  matched  by  the  higher- 
order  linear  systems.  Another  study  by  Wafa  [5]  showed  that  the 
peak  response  executed  by  a  hysteretic  SDF  system  excited  by  an 
earthquake-type  input  is  closely  predicted  by  a  third-order, 
linear,  equivalent  system. 

The  purpose  of  the  present  investigation  is  to  show  that  the 
parameters  of  a  linear,  higher-order,  differential  equation,  ap¬ 
proximately  governing  the  response  of  a  hysteretic  SDF  system, 
can  be  estimated  using  system  identification  techniques.  The 


♦Numbers  in  parentheses  refer  to  a  list  of  references  appended  to 
this  report. 
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parameter  identifications  will  be  based  on  measured  input  (forc¬ 
ing)  and  output  (response)  data.  Parameter  identification  for 
the  linear  system  must  be  attempted  when  the  measured  input  and 
output  signals  are  noisy,  since  some  noise  is  present  in  all 
measured  signals.  Therefore,  system  identification  techniques 
which  can  be  used  successfully  in  the  presence  of  noise  will  be 
demonstrated. 

The  reason  for  pursuing  the  system  identification  study  sum¬ 
marized  in  this  report  is  to  characterize  the  linear  system 
equivalent  to  the  hysteretic  system  in  which  the  damage  has  oc¬ 
curred.  Once  the  linear  system  parameters  are  known,  the  energy 
dissipated  by  this  system  can  be  computed.  It  is  hoped  that  the 
energy  dissipated  by  the  higher-order,  linear  system  is  an  accu¬ 
rate  estimator  of  the  energy  dissipated  by  the  actual  hysteretic 
system.  This  study  addresses  the  question  of  energy  dissipated 
by  these  two  types  of  systems. 

This  report  covers  the  following  subjects:  1)  The  system 
used  to  model  a  hysteretic  SDF  system  is  presented.  This  model 
is  a  higher-order,  linear,  differential  equation.  2)  The  param¬ 
eter  identification  technique  is  given.  Least-squares  parameter 
identification  is  used  in  the  frequency  domain  and  in  the  time 
domain.  3)  Analysis  of  errors  arising  from  the  presence  of  noise 
in  the  measured  inputs  and  outputs  is  presented.  4)  Several  nu¬ 
merical  examples  are  summarized.  5)  Discussion  and  conclusions 
are  given. 

2.0  Model  for  Hysteretic  SDF  System 

The  differential  equation  governing  the  response  of  a  mass 
excited,  potential ly-hysteretic,  single-degree-of-freedom  (SDF) 
structure  is 


where  m  is  the  structure  mass,  f  is  the  forcing  function,  z  is 
the  displacement  response  of  the  mass,  dots  denote  differentia¬ 
tion  with  respect  to  time,  and  u  is  the  restoring  force  of  the 
structure.  In  actual  systems  which  can  be  realistically  modeled 
using  Equation  2.1,  u  can  be  a  complicated  function.  When  the 
system  displacements  are  small,  u  can  reasonably  be  taken  as  a 
linear  function  of  the  displacement  and  velocity;  but  when  the 
displacement  and  velocity  become  large,  u  may  be  a  function  gen¬ 
erating  system  hysteresis  and  whose  characteristics  possibly  vary 
with  time. 

In  analytical  studies,  u  is  modeled  as  a  linear  function  of 
displacement  and  velocity  when  these  quantities  are  small,  and, 
sometimes,  when  they  are  large.  When  an  effort  is  made  to  ac¬ 
count  for  system  hysteresis  accurately,  u  is  taken  as  a  function 
of  time  and  the  response  history.  For  example,  when  a  system  has 
bilinear  hysteretic  stiffness,  its  graph  of  restoring  force  ver¬ 
sus  system  displacement  might  resemble  the  curve  shown  in  Figure 
2.1a.  A  column  made  from  elasto-plastic  material,  with  rectangu¬ 
lar  cross  section,  and  responding  in  bending  will  have  a  restor¬ 
ing  force  versus  displacement  graph  like  the  one  shown  in  Figure 
2.1b. 

We  propose  that  the  response  of  a  system  governed  by  Equa¬ 
tion  2.1  be  modeled  using  the  equations 

mz  +  u  =  f 

(2.2) 

+  z 

j =0 

where  the  cj,  j=0,  ...  Mfl,  are  constants  governing  the  system- 
restoring  force  characteristics ,  u(j)  denotes  the  jth  time  de¬ 
rivative  of  u,  and  M  is  a  constant  denoting  the  order  of  approxi¬ 
mation  of  the  linear  system  to  the  hysteretic  one.  This  model  is 


(j) 

cj  =  CM+1  2 
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DISPLACEMENT 


Figure  2.1 


a.  Spring  restoring  force  versus  displacement  for  bilinear 
nysteretic  structure 


DISPLACEMENT 


Figure  2. 


Spring  restoring  force  versus  displacement  for 
rectangular  column  in  bending 


chosen  to  represent  the  hysteretlc  system  because  it  displays 
hysteresis;  moreover,  it  is  hoped  that  with  the  proper  choice  of 
Cj,  j=0,  ...  M+l,  the  response  of  this  system  will  accurately 
approximate  the  response  of  the  actual  system. 

Note  that  when  M  is  chosen  to  be  0,  the  model  in  Equation 
2.2  becomes 


••  1  *  1 
mz+-iz+-4~z  =  f 

c0  c0 


(2.3) 


This  is  simply  the  linear  differential  equation  which  approxi¬ 
mately  governs  the  response  of  an  SDF  system  whose  response  am¬ 
plitude  is  small.  The  restoring  force  function  for  this  system 
is  u  =  (c  i/c  0)  z  +  (l/c0)  z.  This  system  displays  an  energy  dis¬ 
sipation  characteristic  due  to  viscous  damping.  For  example,  the 
linear  system  whose  parameters  are  listed  in  Figure  2.2a  produces 
the  restoring  force  versus  displacement  graph  shown  in  Figure 
2.2a  when  the  system  is  subjected  to  the  harmonic  force  input 
whose  form  is  given  in  Figure  2.2a.  The  parameters  of  the  model 
in  Equation  2.3  can  be  chosen  so  that  the  model  represents  the 
hysteretic  system  as  well  as  possible,  in  some  sense.  When  the 
response  of  the  actual  system  is  linear  and  damping  is  viscous, 
the  model  of  Equation  2.3  perfectly  represents  the  actual  sys¬ 
tem. 

When  the  constant  M  is  chosen  as  1  in  Equation  2.3,  the 
model  becomes 

mz  +  u  =  f 


CQU  +  CjU  =  c2z  +  z 


(2.4) 


The  restoring  force  is  the  solution  of  the  second  equation.  This 
system  can  display  an  energy-dissipation  characteristic.  For 
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restoring  force  versus  displacement  for  second-order  system 


example,  the  linear  system  whose  parameters  are  listed  In  Figure 
2.2b  produces  the  restoring  force  versus  displacement  graph  shown 
in  Figure  2.2b  when  the  system  is  excited  by  the  harmonic  input 
whose  form  is  given  in  Figure  2.2b.  The  parameters  of  the  system 
in  Equation  2.4  can  be  chosen  so  that  the  model  represents  a  hys- 
teretic  system  as  well  as  possible,  in  some  sense.  When  the 
actual  system  response  is  linear  and  the  system  has  viscous  damp¬ 
ing,  we  anticipate  that  cx  =  0. 

It  is  anticipated  that  as  the  order,  M,  of  the  model  in 
Equation  2.2  is  increased,  the  response  of  a  hysteretic  system 
can  be  matched  with  increasing  accuracy.  However,  for  practical 
reasons  involving  estimation  accuracy  for  system  parameters,  very 
high  order  linear  system  models  cannot  be  used  to  simulate  hys¬ 
teretic  system  behavior. 

Before  considering  the  problem  of  parameter  estimation  for 
the  system  of  Equation  2.2,  note  that  we  anticipate  calculating  a 
set  of  parameters  with  values  in  a  specific  range.  For  example, 
when  M=0  and  Equation  2.3  is  the  model  for  system  response,  we 
anticipate  finding  values  l/c0  >  0  and  cx/c0  >  0  (i.e.,  c0  >  0 
and  cx  >  0).  These  values  guarantee  that  the  model  has  positive 
stiffness  and  damping,  as  we  know  the  real  system  must.  When  M=1 
and  Equation  2.4  models  the  system  response,  we  anticipate  find¬ 
ing  values  c0  >  0,  cx  >  0,  and  c2  >  0.  These  values  guarantee 
that  the  model  response  will  be  stable. 

The  energy  dissipated  by  the  systems  described  in  this  sec¬ 
tion  can  be  computed  using  Equation  2.1  or  2.2.  In  the  terms  of 
Equation  2.1,  the  energy  dissipated  by  a  system  is 


Figure  2.2b.  Total  restoring  force  versus  displacement  for  third-order  system 


where  z(0)  is  the  system  displacement  at  time  0,  z(T)  is  the  sys¬ 
tem  displacement  at  time  T,  and  T  is  the  time  through  which  the 
energy  computation  is  executed.  This  formula  provides  the  area 
enclosed  in  the  hysteresis  loops  generated  by  the  system  re¬ 
sponse.  For  example,  Equation  2.5  could  be  used  to  compute  the 
area  enclosed  by  the  hysteresis  loop  of  Figure  2.1b.  This  equa¬ 
tion  can  be  transformed  into  terms  more  convenient  for  computa¬ 
tion.  Note  that  the  integral  of  Equation  2.5  is  written  in  terms 
of  the  displacement  variable  z.  The  variable  of  integration  can 
be  transformed  to  a  time  variable  yielding  the  following  expres¬ 
sion. 


-D 


J  u  z  dt 


(2.5a) 


Finally,  Equation  2.1  can  be  used  above  to  obtain 
(f  -  mz)  z  dt 


■r 


(2.5b) 


When  the  input  and  response  for  an  SDF  system  are  measured, 
Equation  2.5b  can  be  used  to  directly  compute  the  energy  dissi¬ 
pated  by  a  system.  When  the  response  is  computed,  for  example  by 
solving  Equations  2.4,  then  the  input  and  system  parameters  are 
used  to  find  z  and  its  derivatives.  The  input,  f,  and  computed 
response,  z,  are  used  in  Equation  2.5b  to  determine  the  energy 
dissipated. 

Numerical  examples  where  we  compute  the  parameters  of  linear 
systems  equivalent  to  bilinear,  hysteretic  systems  are  presented 
in  Section  5.  In  these  examples,  the  energy  dissipated  by  each 
system  is  computed  and  the  results  are  compared. 


3.0  Identification  of  Parameters 

We  now  consider  the  problem  of  parameter  identification  for 
the  linear  model  for  hysteretic  single-degree-of-freedom  (SDF) 
systems  introduced  in  the  last  section.  Both  the  time  and  fre¬ 
quency  domain  approaches  to  the  problem  of  parameter  identifica¬ 
tion  are  outlined.  We  consider  the  parameter  identification  both 
when  noise  is  and  is  not  present  in  the  measured  data. 

Measured  input  and  output  data  from  a  structural  system  are 
used  to  estimate  the  system  parameters.  It  is  assumed  that  these 
data  are  available.  To  begin  with,  we  assume  that  no  noise  is 
present  in  the  measured  data.  Later,  modifications  are  made  to 
permit  the  inclusion  of  noise.  The  input,  or  forcing  function, 
data  are  available  in  discrete  form  and  are  denoted  f4,  i=0, 

...  n-1.  The  output,  or  system  response,  data  are  available  in 
discrete  form  and  are  denoted  zA,  t=0,  ...  n-1.  It  is  assumed 
that  the  input  and  output  signals  have  been  sampled  at  times  t^ 

=  fcAt,  *=0,  ...  n-1.  Note  that  measured  response  data  are  as¬ 
sumed  to  be  given  as  accelerations;  this  is  realistic  since 
structural  response  acceleration  is  often  the  quantity  measured 
during  an  experiment. 

The  approaches  to  parameter  identification  developed  in  this 
section  do  not  require  that  the  measured  data  be  indexed  by 
time.  Any  sampled  sets  f*  and  zA,  z=0,  ...  n-1  can  be  used. 
However,  since  the  present  applications  come  from  the  field  of 
structural  dynamics,  we  will  refer  to  signals  sampled  in  time. 

The  model  established  in  Equation  2.2  is  repeated  here;  it 
is 

mz  +  u  =  f 

(3.!) 

Cj  U<J  *  CM+1  *  *  2 
j-0 
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When  the  first  of  these  equations  is  solved  for  u  and  the  result¬ 
ing  expression  is  used  in  the  second  equation,  we  obtain 


(3.2) 

where  f^)  is  the  jth  derivative  of  f,  and  z^+2^  is  the  (j+2)th 
derivative  of  z.  (This  particular  arrangement  is  chosen  since 
c0  should  never  be  0  for  the  systems  under  consideration.)  The 
notation  in  this  equation  can  be  simplified  by  taking 


M. 


(3.3) 


The  equation  governing  motion  of  the  system  becomes 
M 

aQz  +  a^z  +2^aj+i  (_f^  +  mz^J+2^)  =  f  -  mz. 
j=l 

(3.4) 

This  equation  can  be  used,  along  with  the  measured  input  and  re¬ 
sponse  data,  to  identify  the  parameters  a j ,  j=0,  ...  M+l.  When 
these  are  known.  Equation  3.3  can  be  used  to  find  the  cj,  j=0, 
...  M*l. 

3.1  Time  Domain  Approach  I 

Let  z^,  z^,  and  z  t=0,  ...  n-1,  be,  respectively,  the 
response  displacement,  velocity  and  jth  derivative  of  the  dis¬ 
placement  at  time,  t^,  z=0,  ...  n-1.  Let  f i  and  f ^ ,  t=0, 

...  n-1,  be  the  force  at  time  tt,  1=0,  ...  n-1.  Assume  that 
all  required  integrals  and  derivatives  of  and  z4  can  be 
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obtained,  exactly.  Then  Equation  3.4  can  be  written  at  time  tA 


ao2i  +  a 


1*1  ”Zt(J*2))  *  \  -  «t.  »■».  ...  n-l 


(3.6) 


The  notation  in  this  equation  can  be  simplified  by  defining 
{z^}  ■  {zl  z l  -ft  +  mzt  ...  +  mz^M+2^)  ,  4=0,  ...  n-l 


{a }  *  Ug  ai  •••  aM+i) 


(3.6a) 


(3.6b) 


where  the  T  superscript  refers  to  matrix  transposition.  Use  of 
these  expressions  in  Equation  3.5  yields  the  relation 


(Z,JT  (al  =  4*0,  ...  n-l 


(3.7) 


This  is  the  equation  governing  the  linear  system  response  at  time 
tf 

The  notation  can  be  further  simplified  by  defining 


WT 


f0  '  mz0 


[zf ]  -  {zi]T  ,  {fz}=  )  fi  ' 


(3.8) 


lzn-ltT 


fn-r  mzn-l 


Using  these  expressions,  the  sequence  of  Equation  3.7,  for  4=0, 
...  n-l,  can  be  written 


[Zf]  {a}  -  {fz} 


(3.9) 


This  equation  governs  the  linear  system  response  at  all  times. 

When  1)  the  mechanical  system  from  which  the  data  were  mea¬ 
sured  is  truly  linear,  2)  there  is  no  noise  in  the  data  measure¬ 
ments,  and  3)  the  derivatives  and  integrals  of  and  z^,  l= 0, 

...  n-1,  are  known  exactly,  Equation  3.9  can  be  satisfied  exactly 
by  the  measured  data  and  a  set  of  coeff icients ,  aj,  j=0,  ... 

M+l.  Otherwise,  the  left-hand  side  of  Equation  3.9  will  not 
exactly  equal  the  right-hand  side  at  every  time,  t£.  We  define 
the  error  vector  as  follows. 

[£1  =  ( eQ  ei  •••  £n_i)  (3.10) 

-  ^f]  [a]  -  [f2] 

The  elements  e£,  1=0,  ...  n-1,  quantify  the  data  nonlinearity, 
the  measurement  noise,  and  the  inaccuracy  of  the  derivatives  and 
integrals  of  z"A  and  f^,  as  a  function  of  the  system  parameters, 
{a}. 

An  overall  measure  of  the  lack  of  fit  of  the  measured  data 
to  the  model  in  Equation  3.9  is  established  in  the  parameter 

s2-  |«|T{.|-  ((alT[ZfiT-  {<■/)  ([*,](•!-  (f2l) 

(3.11) 

This  is  referred  to  as  the  squared  error  between  the  measured 
data  and  the  system  model.  This  error  can  be  minimized  through 
the  proper  choice  of  the  parameter  vector  {a}.  The  vector  {a} 
which  minimizes  e2  is 

{*}  «  C  [Zf  ]T[Zf  ])  1  [Zf]T|fz}  (3.12) 
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(See,  for  example.  Reference  6.)  The  parameter  vector  chosen 
above  is  the  best  estimate  of  the  system  parameters  in  a 
least  square  sense. 

If  the  quantity,  e2/n,  is  small  relative  to  the  signal  dif¬ 
ference,  f  -  mz^,  l=0,  ...  n-1,  then  the  model  of  Equation  3.4 
with  parameters  obtained  using  Equation  3.12  accurately  repre¬ 
sents  the  measured  system.  If  this  is  not  true,  then  the  model 
of  Equation  3.4  is  inadequate. 

e2  will  equal  zero  if  the  measured  system  data  are  noise 
free,  the  system  is  linear,  and  the  computed  derivatives  and  in¬ 
tegrals  of  the  data  are  exact.  Failure  to  meet  any  one  or  more 
of  these  requirements  will  cause  e2  to  be  nonzero.  In  practice, 
the  parameter  identification  procedure  outlined  above  can  only  be 
used  effectively  when  there  is  little  or  no  noise  in  the  measure¬ 
ments,  the  system  is  nearly  linear,  and  the  integrals  and  deriva¬ 
tives  of  z  and  f  can  be  obtained  accurately.  The  method  is  par¬ 
ticularly  effective  when  the  parameter  M  is  set  to  0.  However, 
when  noise  is  present  and  M=l,  the  procedure  loses  accuracy  since 
derivatives  of  z  and  f  are  required.  Particularly,  the  functions 
z  and  f  must  be  estimated.  If  the  measured  raw  data,  iz  and 
f^,  t=0,  ...  n-1,  are  used  directly  to  obtain  these  estimates 
through  numerical  differentiation,  the  estimates  of  z  and  f  may 
be  very  poor  due  to  amplification  of  the  effects  of  noise  and  in¬ 
accuracy  of  numerical  differentiation.  In  view  of  this,  alter¬ 
nate  procedures  for  parameter  identification  in  the  presence  of 
noise  must  be  established. 

It  is  possible  to  retain  a  time  domain  approach  in  the  pa¬ 
rameter  estimation  problem.  For  example,  the  signals  iz  and 
fjj,  t=0,  ...  n-1  can  be  smoothed  prior  to  the  taking  of  deriv- 
ctives.  Smoothing  can  be  accomplished  through  the  use  of  digital 
filters;  or  the  governing  Equation  3.4  can  be  integrated  M  times 
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before  applying  the  identification  procedure.  Then  no  numerical 
derivatives  are  required.  Alternately,  a  frequency  domain  ap¬ 
proach  can  be  used.  This  eliminates  the  necessity  to  directly 
differentiate  measured  data.  In  the  present  study,  we  pursue  the 
approach  using  the  integrated  equations  of  motion  and  the  fre¬ 
quency  domain  approach.  First,  we  consider  the  time  domain  ap¬ 
proach  using  the  integrated  equations  of  motion. 

3.2  Time  Domain  Approach  II 


Consider  again  Equation  3.4.  That  equation  establishes  the 
assumed  form  for  the  governing  equation  of  motion  of  a  hysteretic 
SDF  system.  The  direct  use  of  that  equation  in  a  parameter  iden¬ 
tification  scheme  requires  the  calculation  of  z(M+2)(t),  the 
Mth  derivative  of  z(t).  The  numerical  computation  of  derivatives 
is  an  operation  which  is  generally  sensitive  to  the  presence  of 
noise  in  measured  signals;  therefore,  we  employ  the  following 
procedure.  Integrate  both  sides  of  Equation  3.4  M  times.  We 


obtain 


a0*M  +  al*M-I  +  E  Vl  (-gM-j  + 

j  =  l 


=  gM  -  *yM_2 


(3.13) 


where  y_2(t)  =  z(t) 


=  f  yj-i <T)  dt 

J  n 


j = - 1 ,  ...  M 


gn(t)  ■  f(t) 


Sj(t)  =r  gj. 

•'n 


(*)  dl 


j=l,  ...  M. 


It  is  assumed  that  the  initial  values  of  z(t)  and  f ( t )  and  their 
derivatives  are  zero.  Equation  3.13  involves  ’z(t)  and  its  first 
M+2  integrals,  and  f(t)  and  its  first  M  integrals.  When  ’z(t)  and 
f(t)  are  measured,  these  can  be  calculated  numerically. 

A  parameter  identification  can  now  be  executed  using  Equa¬ 
tion  3.13.  As  before,  it  is  assumed  that  z(t)  and  f(t)  have  been 
measured  at  discrete  times,  t4  =  tAt,  1=0,  ...  n-1.  The  inte¬ 
grals  of  z(t)  and  f(t)  can  be  directly  computed  at  the  times 
t^.  At  time  t*tA,  y j (t )  is  denoted  y^  v  g j (t )  is  denoted  g^ 
and  Equation  3.13  becomes 


M 

a0yM,z  +  alyM-l,z  +  £aj  +  l(_gM-j,*  +  =  9M,t  "  myM-2,t 

j  =  l 

A=0,  ...  n-1 
(3.14) 


This  equation  governs  the  SDF  system  at  each  time  t4.  The  no- 


Pfl 


for  the 

collection  of  equations  described  i 

n  3.14 

can  be 

fied  by 

defining 

the  matrices 

(a0 

•  •  •  3 

)T 

M+l ' 

C 

!.15a) 

yM,0 

yM-l ,0 

^*gM-l  ,0 

<-go,o 

*"y-2,0> 

yM,l 

• 

yM-l ,1 

• 

(_gM-l  ,1  'W,yM-  3 , 1  * 

• 

(-go,i 

^-2,1) 

• 

• 

• 

yM,n-l 

• 

• 

yM-l,n- 

• 

• 

1  ^■gM-i,n-l+myM-3,n-l^  *** 

(-g0,n- 

• 

• 

l+my-2,n-l'' 

(3.15b) 

“  ^gM,0'  myM-2,0^  (gM,r  myM-2,P  ^M.n-l"  myM-2,n-l^ 

(3.15c) 
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In  terms  of  these  matrices.  Equation  3.14  can  be  expressed 


[Zf]  (a!  =  {fzl  (3.16) 

This  equation  has  the  same  form  as  that  used  in  the  parameter 
identification  problem  for  Equation  3.4,  but  the  matrices  [Zf] 
and  {f2 }  have  been  defined  differently.  When  1)  a  physical 
system  is  exactly  governed  by  the  (M+2)th  order  Equation  3.1,  2) 
there  is  no  noise  in  the  measurements  of  z(t)  and  f(t),  and  3) 
the  integrals  of  z(t)  and  f ( t )  are  computed  exactly.  Equation 
3.16  can  be  satisfied  exactly  by  the  measured  data  and  a  collec¬ 
tion  of  coefficients,  aj,  j=0,  ...  M+l.  However,  when  these 
requirements  are  not  satisfied,  Equation  3.16  will  not  hold, 
exactly.  Then  the  equation 

M  =  [Zf]  {*}  -  [f2]  (3.17) 

will  hold,  where  { e }  is  an  error  vector  which  reflects  the  degree 
to  which  the  above  requirements  are  not  satisfied.  As  before, 
the  sum  of  the  squares  of  the  components  in  { e }  is  minimized  by 
choosing 

(a  I  •  <[Zf]T[Zfr!  [Zf]T(f2l  (3.18) 

When  the  elements  in  {e}  are  relatively  small,  Equation  3.1 
provides  a  satisfactory  model  for  the  system  under  consideration 
and  the  vector  |a }  includes  the  parameters  of  that  model. 

The  advantage  to  the  present  approach  to  parameter  identifi¬ 
cation  is  that  the  derivatives  of  z(t)  and  f(t)  need  not  be  com¬ 
puted  via  numerical  differencing.  This  is  a  practical  improve¬ 
ment  since  all  measured  signals  contain  at  least  some  noise. 
Beyond  this,  the  operation  of  numerical  differencing  introduces 


so  much  Inaccuracy  Into  a  computation  that  parameter  identifica¬ 
tion  is  difficult,  even  when  the  signals  contain  no  noise.  The 
use  of  Equation  3.13  replaces  the  requirement  for  differentiation 
with  an  integration  requirement.  When  the  measured  signals  are 
noisy,  their  integrals  are  influenced  by  the  presence  of  noise; 
but  the  noise  influence  is  less  than  in  the  case  of  differentia¬ 
tion.  In  Section  4,  the  effect  of  noise  on  parameter  identifica¬ 
tion  will  be  discussed. 

3.3  Frequency  Domain  Approach 

The  final  approach  we  will  consider  for  the  identification 
of  the  parameters  of  Equation  3.4  is  a  frequency  domain  ap¬ 
proach.  Consider  Equation  3.4.  Both  sides  can  be  Fourier  trans¬ 
formed  to  obtain 


M 


a q  Z(w)  +  iwZ(w)  +  Z>I..  (-(iw)^F(w)  +  m(iw)^+^Z(w) ) 

j=l 


=  F(w)  +  mw^Z(w) 


where  Z( 


r 

—  a 


z(t)  e'iwtdt 


—  <  w  <  -  (3.19) 


-«  <  w  <  <® 


fm  .  r 

•'-at 


f(t)  e'1wtdt 


-«  <  w  <  • 


are  the  Fourier  transforms  of  z ( t )  and  f(t),  respectively.  This 
equation  can  be  rearranged  by  grouping  the  coefficients  of  Z(w) 
and  F(w),  then  combining  Z(w)  and  F(w)  on  one  side  of  the  equa¬ 
tion. 


-mw 


2 


+  a, 


+  aj(iw)  +  m^aj+1(iw)j+2 

j=l 

- R - 

1  +Zvi(iw)j 

j=l 


-*>  <  w  <  » 


(3.20) 


Finally,  each  side  of  the  equation  can  be  multiplied  by  its  com¬ 
plex  conjugate  to  yield  the  modulus  squared. 


1*1 

|-mw2  +  aQ  +  a  x  (i  w)  +  m  ^TV+1  (i  w)J+2[ 
_ _ ill _ 

1  +£vi«">J|2 

j-1 


F(w) 


Z(w) 


2* 


<  w  < 


(3.21) 


When  this  expression  is  evaluated  at  w  =  w^,  the  result  is 


=  ’ *  * 


(3.22) 


where  | | 2/|zk| 2  has  been  set  equal  to  Q^.  It  is  assumed  that 
the  estimates  for  the  Fourier  transforms  of  z(t)  and  f(t)  can  be 
obtained  by  using  the  discrete  Fourier  transform  (D FT)  of  zA 


21 


C  |C\J 


and  fA,  t=0,  ...  n-1.  Then  these  estimates  will  be  available 
at  the  frequencies  W|<  =  2irk/nat,  k=0,  ...  n/2.  The  DFT  can  be 
efficiently  obtained  using  the  fast  Fourier  transform  (FFT)  pro¬ 
cedure.  (See,  for  example,  Reference  7).  The  DFT  of  zit  z=0, 

...  n-1,  is  assumed  to  be  -wk2  Zk ,  k=0,  ...  n/2.  Zk  and  Fk ,  k=0, 
...  n/2,  represent  the  DFTs  of  z(t)  and  f(t)  at  frequencies  W|<, 
k=0,  ...  n/2. 

When  1)  the  system  under  consideration  is  linear  and  or 
order  M+2,  2)  all  measurements  are  noise  free,  and  3)  the  DFTs  of 
z l  and  f^,  ...  n-1  used  in  Equation  3.22  are  exact,  Equation 

3.22  can  be  satisfied  exactly  by  the  measured  data  and  a  set  of 
coefficients,  aj,  j=0,  ...  M+l.  When  these  requirements  are 
not  satisfied.  Equation  3.22  cannot  be  precisely  satisfied.  In 
this  latter  case,  we  consider  the  equation 


-"■"k2  * 


ao  *  ai«"k>  *  »£aj+i««k>J+2 

_ _ _ iil _ 

R ~2 

I1  ^aj*l(1uk,J  I 


=  Qk  £k»  k=0,.. . 


(3.23) 


On  the  right-hand  side  is  an  error  term.  This  term  measures 
the  degree  to  which  the  model  of  Equation  3.4  represents  or  fails 
to  represent  the  measured  data,  when  the  measured  data  are  used 
to  compute  F|<  and  Z|<,  k=0,  ...  n/2.  We  wish  to  find  the  co¬ 
efficients  aj,  j=0,  ...  M+l,  which  minimize  the  e^,  k=0,  ... 
n/2.  This  can  be  accomplished  using  a  least-squares  criterion 
where  a  sum  £e)<  2  is  minimized.  This  can  be  done  numerically  by 
searching  through  the  space  of  (M+2) -tuples  created  by  the  coor¬ 
dinates  ag,  ...  For  example,  a  gradient  search  technique 

can  be  used  to  find  the  minimum.  When  M  is  chosen  as  0  or  1  and 
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when  the  system  to  be  analyzed  has  a  response  which  is  only  mod¬ 
erately  nonlinear,  then  a  simpler  technique  for  choosing  the  best 
aj,  j=0,  ...  Mfl,  is  available.  Consider  first  the  case  where 
M=0.  When  M  is  chosen  as  zero.  Equation  3.4  becomes  the  equation 
of  a  linear,  second-order,  SDF  system.  The  equation  to  be  used 
in  the  parameter  identification.  Equation  3.23,  can  be  written 

2  2 

J-mwk  +  aQ  +  a1  iwk|  =  Qk  +  ek>  k=0,  ...  n/2  (3.24) 

The  left  side  of  the  expression  can  be  expanded  to  yield 

(aQ  -  mwk2)2  +  (a!w|<)2  =  Qk  +  ek.  k=0»  ...  n/2  (3.25) 


The  functional  expression  on  the  left-hand  side  assumes  its  mini¬ 
mum  when  w  takes  the  value 


w 


1 

m 


<ao  -  sr> 


1/2 


(3.26) 


Note  from  Equation  3.4  that  a0  is  the  equivalent  stiffness  and 
a!  is  the  equivalent  damping  for  a  system  whose  Equation  3.4 
model  uses  M=0.  In  view  of  this,  whenever  the  equivalent  damping 
factor  is  much  less  than  1  (say  less  than  0.2),  aj2/2m  is  small 
compared  to  a0.  (This  will  be  true  in  most  cases.)  In  this 
case,  a0  can  be  estimated  by  noting  the  frequency  where  Qk  + 
e](,  k=0,  ...  n/2  is  a  minimum.  Denote  this  frequency  wm  and 
use  it  in  Equation  3.25  to  show  that 


a0  = 


(3.27) 


To  obtain  Improved  accuracy,  the  frequency  where  Qk  +  ek,  k=0, 
...  n/2  is  a  minimum  can  be  estimated  by  using  the  least-squares 
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approach  to  fit  a  parabola  to  the  values  Qk  +  ek  near  the  char¬ 
acteristic  frequency,  .  The  parameters  of  this  parabola  can 
be  used  to  determine  the  frequency  where  Qk  +  ek  is  a  minimum. 
Estimation  of  the  constant  a0  by  any  means  must  assume  that  no 
quality  in  ek*  k=0,  ...  n/2,  will  tend  to  move  the  minimum 
value  of  Qk,  k=0,  ...  n/2,  in  a  systematic  way.  The  general 
character  of  ,  k-0,  ...  n/2,  will  be  discussed  and  modeled  in 
Section  4. 

The  motivation  for  using  the  above  analysis  to  estimate  a0 
is  clarified  by  the  following  considerations.  The  characteris¬ 
tics  of  an  SDF  system  are  reflected  in  the  system  behavior  at  all 
frequencies  but  are  best  interpreted  from  its  behavior  near  the 
characteristic  frequency.  (For  strictly  linear  response,  we 
would  refer  to  the  natural  frequency,  but  the  response  considered 
here  is  not  always  linear.)  The  frequency  band  including  the 
characteristic  frequency  is  that  band  where  the  system  amplifies 
the  response,  and  in  this  band  the  values  of  Qk ,  k=0,  ...  n/2, 
form  a  trougn.  Thfs  is  the  reason  why  the  minimum  value  of  Qk 
+  ek»  k=0,  •••  n/2,  is  considered,  above. 

The  estimated  value  of  a0  can  now  be  used  in  Equation  3.25 
to  obtain 

m2(wm2  -  wk2)2  +  (a^wk)2  =  Qk  +  ek,  k=0,  ...  n/2  (3.28) 

Now,  in  order  to  estimate  aj,  we  can  employ  a  least-squares  cri¬ 
terion.  Specifically,  we  solve  Equation  3.28  for  ck»  square 
this  expression,  and  sum  over  all  or  part  of  the  values  of  k. 

The  resulting  expression  is  denoted  e2.  This  quantity  is  mini¬ 
mized  with  respect  to  a^  by  differentiating  e2  with  respect  to 
a!2,  equating  the  result  to  zero,  and  solving  for  a^  The  opera¬ 
tions  listed  above  yield  the  result 


(3.29) 


al  =  Wk  ^  “  m2^wm2  " 

k  k 

The  sums  in  Equation  3.29  are  listed  simply  as  sums  over  k.  In 
the  identification  process,  these  sums  should  be  taken  over  those 
frequencies  near  the  characteristic  frequency  of  the  SDF  system. 
These  are  the  frequencies  where  the  ratio  |  | 2/ | | 2  best  char¬ 

acterizes  the  properties  of  the  SDF  system.  For  example,  the 
values  of  k  used  in  the  analysis  might  be  chosen  as  those  corre¬ 
sponding  to  the  frequencies,  wk,  falling  in  the  interval  which 
is  twice  the  bandwidth  of  the  SDF  system  under  consideration.  Of 
course,  this  band  of  frequencies  will  not  be  precisely  known 
since  the  system  parameters  have  not  been  identified.  However, 
it  can  be  estimated  as  the  frequency  band  (w.  ,  w.  ).  w.  is  the 

,  I?  I  y  M  k2  k1 

lower  frequency  where  |  Ffc  \  /\\  is  about  five  times  as  great 
as  the  minimum  value  of  |Fk|2/|Zk  2.  wfe  is  the  higher  frequency 
where  |Fk  |2/|zk  |2  is  about  5  times  as  great  as  the  minimum 
value  of  1%  j2/|Z|<  |2. 

The  constants  a0  and  a^  expressed  in  Equations  3.27  and 
3.29  provide  the  parameters  for  a  second-order  model  for  an  SDF 
system.  This  model  is  the  best  available  in  a  least-squares 
sense.  Note  again  that  the  means  used  to  obtain  a0  assumes  rela¬ 
tively  low  damping. 

Some  may  consider  it  undesirable  or  inaccurate  to  model  a 
hysteretic  SDF  system  using  a  second-order,  linear  ordinary  dif¬ 
ferential  equation,  as  we  do  with  Equation  3.4  when  M=0.  With 
the  hope  of  obtaining  improved  accuracy,  we  can  model  the  SDF 
system  under  consideration  using  Equation  3.4  with  M=l.  Note 
first,  though,  that  1)  improved  accuracy  in  every  sense  may  not 
be  obtained  using  M=l.  This  is,  using  M=1  in  Equation  3.4  only 
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guarantees  that  the  parameters  aj,  j=0,l,2  will  yield  a  more 
accurate  match  between  the  motion  response  of  the  model  and  the 
measured  response;  the  energy  dissipation  characteristics  of  the 
model  system  may  not  be  more  accurate.  2)  The  linear  viscous 
damping  term  used  in  most  linear  system  analyses  does  not  pre¬ 
cisely  model  real  system  behavior,  even  when  the  real  system  has 
small  displacements.  Yet,  linear  viscous  damping  is  freely  used 
in  analyses. 

When  M  is  chosen  as  1  in  Equation  3.4,  the  formula  to  be 
used  in  the  parameter  identification  is  Equation  3.23  with  M=l. 
This  is 


-mwk2  +  aQ  +  a1(iwk)  +  ma2(iwk) 
- - - 

1  +  a2(iwk) 


3.2 


Qk  £k»  k=0,  • 


..  n/2. 

(3.30) 


The  numerator  and  denominator  can  be  expanded  on  the  left  side  to 
obtain 


2  2  2  2  2 
(an  -  mwu  )  +  wt  (a1  -  mwk  a2) 


1  +  (a2wk)' 


~  Qk  +  ek ,  k=0, 


n/2. 


(3.31) 


When  a2  is  small  compared  to  the  characteristic  frequency,  wm 
of  the  SDF  system  (say  a2  <.  0.3  v^) ,  the  (a^)2  term  in  the 
denominator  can  be  neglected.  This  is  usually  the  case  when  non¬ 
linear  deformations  are  not  too  large.  When  a!  and  a2  are  small 
compared  to  a0  (which  is  usually  the  case)  the  minimum  of  the 
left-hand  side  occurs  near  the  frequency 
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(3.32) 


We  find  this  frequency  by  scanning  Qk  +  ek»  k*0,  . ..  n/2 ,  or 
by  using  the  curve  fitting  technique  discussed  above.  This  fre¬ 
quency  is  denoted  %.  In  terms  of  a0  can  be  written 


a 


0 


(3.33) 


This  expression  and  the  previous  assumptions  can  be  used  in  Equa 
tion  3.31  to  obtain 

m  (w^  “  (a^  -  mwk  a2)  =  Qk  +  £k»  k=0,  ...  n/2 

(3.34) 

When  mwm2a2  is  relatively  small  compared  to  ax,  the  second  term 
on  the  left  side  of  Equation  3.34  can  be  expanded,  and  the 
m2Wk4a22  term  can  be  neglected.  This  is  possible  when  the 
amount  of  yielding  in  the  system  to  be  modeled  is  not  too  great. 
The  result  of  this  operation  is 

2 

m2(wm2  -  wk2)  +  wk2(a12  -  2mwk2a1a2)  =  Qk  +  ek,  k=0,  ...  n/2 

(3.35) 


Finally,  we  can  set 
2 

al  =  bl  and  ala2  =  *32  (3.36) 

to  obtain 

2 

wJ(2(b1  -  2mwk2  b2)  -  (Qk  -  m2(wm2  -  wfc2)  )  =  ek,  k=0,  ...  n/2 

(3.37) 
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The  system  parameters,  ai  and  a2>  can  now  be  identified  as  those 
which  minimize  the  sum  of  the  squares  of  the  terms  in  Equa¬ 
tion  3.37. 

A  notational  simplification  can  be  achieved  for  the  sequence 
of  Equations  3.37  by  defining 


{b}  =  (bj  b2)T 


w 


^J+l 


(3.38a) 


(3.38b) 


(3.38c) 


(3.38d) 


Note  that  here,  as  in  the  identification  of  the  parameters  of  the 
second-order  system,  only  a  portion  of  the  Equations  3.37  will  be 


used  in  the  parameter  identification.  In  particular,  the  identi¬ 
fication  is  performed  in  the  frequency  range  (w.  ,  w,  ).  This 

K 1  *2 

frequency  range  can  be  chosen  as  before,  using  the  guidelines 
presented  following  Equation  3.29.  In  terms  of  these  matrices. 
Equations  3.37  can  be  expressed 

[XJ  {b}  -  {X2}  =  {e}  (3.39) 

We  wish  to  find  the  vector  (b|  which  minimizes  the  constant  e2  = 
WTW-  it  can  be  shown  (see,  for  example,  Reference  6)  that 
this  is 


-1 


{b  }  =  ( [X1]T[X1])  [X1]T{X2} 


(3.40) 


When  bx  and  b2  have  been  computed,  a!  and  a2  can  be  found  from 


al 


a2  =  b2/a1 


(3.41) 


The  constants  ao,  als  and  a2,  defined  in  Equations  3.33  and 
3.41,  are  the  parameters  of  the  third-order  linear  system.  Equa¬ 
tion  3.4,  whose  response  best  matches  the  hysteretic  response 
z(t),  defined  earlier  in  this  section. 

Some  numerical  examples  are  presented  in  Section  5.  These 
demonstrate  the  use  of  all  the  parameter  identification  tech¬ 
niques  outlined  above.  Moreover,  the  examples  compare  the  energy 
dissipated  by  the  equivalent  systems  of  Equation  3.4  to  the  en¬ 
ergy  dissipated  in  some  actual  nonlinear  systems. 

4.0  Error  Analysis 

In  the  previous  section,  two  techniques  used  to  establish 
the  model  parameters  for  hysteretic  systems  were  presented.  Both 
time  domain  and  frequency  domain  approaches  for  parameter  estima¬ 
tion  were  outlined.  At  various  points  in  Section  3,  allusions  to 
the  potential  for  presence  of  noise  in  the  measured  input  and 
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output  signals  were  made.  The  approaches  chosen  for  use  in  pa¬ 
rameter  identification  were  chosen  for  their  capacity  to  accu¬ 
rately  estimate  system  parameters  in  the  presence  of  noise. 

In  this  section  we  show  how  an  approximate  analysis  of  the 
accuracy  of  the  estimated  parameters  can  be  executed,  when  the 
parameters  are  estimated  using  the  frequency  domain  approach. 

The  presence  of  noise  in  the  measured  inputs  and  outputs  affects 
this  accuracy. 

In  this  analysis,  the  following  assumptions  are  made.  1) 
There  exists  in  nature  a  true  input,  f(t),  and  a  true  structural 
response,  z(t).  These  functions  possess  exact  derivatives  and 
integrals,  and  if  these  were  known,  then  any  of  the  parameter 
identification  techniques  given  in  Section  3  could  be  used  to  es¬ 
tablish  the  parameters  in  the  model.  Equation  3.4,  exactly.  2) 
The  measured  structural  input  and  response  may  contain  a  noise 
component,  in  practice.  When  noise  is  present  in  a  measured  in¬ 
put  or  response,  it  appears  as  an  additive  term.  When  noise  is 
present  in  the  measured  input  or  response,  it  is  zero-mean, 
normally-distributed,  band-limited  white  noise.  Noise  signals 
added  to  the  input  and  output  are  independent.  When  the  struc¬ 
tural  input  and  output  signals  are  discretized  at  a  time  interval 
of  At,  the  cutoff  frequency  for  the  noise  signals  is  l/2At 
hertz. 

In  the  error  analysis  to  follow,  we  will  explicitly  consider 
the  effects  of  measurement  noise  on  parameter  identification. 

The  effects  of  nonlinearity  of  the  measured  system  response,  and 
accuracy  of  the  Fourier  transform  will  not  be  considered. 

Recall  from  Section  3  that  when  the  frequency  domain  ap¬ 
proach  is  used  to  estimate  the  parameters  aj,  j  =  0,  ...  M+l, 
these  parameters  may  be  estimated  using  a  two-step  procedure.  In 
the  first  step,  a0  is  estimated,  and  then  the  remaining  param¬ 
eters,  aj ,  j=l,  ...  M+l,  are  estimated.  a0  is  estimated  by 
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noting  the  minimum  value  in  +  e^,  k=0,  ... 
is  defined 


n/2,  where  Qk 


k=0,  ...  n/2 


(4.1) 


We  now  consider  the  effect  that  the  presence  of  noise  on  the  in¬ 
put  and  output  have  on  Q^. 

In  Action  j  it  was  stated  that  the  input  forcing  function 
f t=0,  ...  n-1  and  the  response  acceleration,  z^,  4=0,  ...  n-1 
are  assumed  known.  When  the  measured  forcing  function  and  re¬ 
sponse  include  noise,  these  functions  become  f  +  n^,  4=0,  ... 
n-1,  and  z ^  +  s^,  4=0,  ...  n-1,  respectively,  and  s^,  4=0, 

...  n-1,  are  the  input  and  response  noise  random  processes. 

These  are,  as  stated  above,  zero-mean,  normally-distributed, 
white  noise,  random  processes.  These  are  discrete  time  random 
processes.  The  white  noise  assumption  implies  that  all  random 
variable  n  ,  n  are  independent  for  l*r;  s  and  s  are  indepen- 

Xf  I  Xr  I 

dent  for  4*r.  The  random  processes  nj,  4=0,  ...  n-1,  is 
independent  of  s^,  4=0,  ...  n-1.  The  variances  of  the  random 
processes  n^  and  s^,  4=0,  ...  n-1  will  be  denoted  an2  and  a$2, 
respectively. 

Since  n^and  s^,  4=0,  ...  n-1  are  normally-distributed,  ran¬ 
dom  processes,  they  are  fully  characterized  by  their  first-  and 
second-order  moments.  The  mathematical  expressions  which  sum¬ 
marize  the  assumptions  listed  above  follow.  The  zero  mean  as¬ 
sumption  is  expressed 

E[nJ  *  e[s4]  =0,  4=0,  ...  n-1  (4.2) 

where  E [ •  ]  refers  to  the  operation  of  expectation.  The  autocor¬ 
relation  functions  of  the  random  processes  n^  and  s^,  4=0,  ... 
n-1  reflect  the  white  noise  assumption.  These  are 
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ti- -a. 


i,  r=0,  ... 


n-1 


(4.3a) 


r 


£(",%]  ■  °n  hr  • 

E[s,  h>  '  «s2  V  •  *•  r-°. 


n-1 


(4.3b) 


where  6^  is  the  Kronecker  delta. 

To  execute  a  frequency  domain  analysis,  we  must  find  the 
discrete  Fourier  transforms  (DFT)  of  f^  +  and  z  +  s^,  i= 0, 
...  n-1.  The  DFTs  of  f^  and  z^,  i-O,  ...  n-1,  have  been  ex¬ 
pressed  F^  and  -w2Zk,  k=0,  ...  n/2,  respectively.  The  DFTs  of 
n^  and  s^,  t=0,  ...  n-1,  are  defined 


Nk  =  e_l2irik/n»  k=0»  •••  n/2 


-i  2irtk/n 


k=0. 


n/2 


(4.4a) 


(4.4b) 


The  sequences  N|<  and  S^,  k=0,  ...  n/2,  form  normally  distrib¬ 
uted  random  processes  since  each  of  the  values  and  S^,  is 
a  linear  function  of  a  sequence  of  normally-distributed  random 
variables.  In  view  of  this,  and  S^,  k=0,  ...  n/2,  can  be 
fully  characterized  in  terms  of  their  first-  and  second-order 
moments.  Those  are  listed  below. 

Each  of  the  random  variables,  \ ,  S^,  is  a  complex  quan¬ 
tity,  and  possesses  a  real  and  an  imaginary  part.  In  terms  of 
these,  N|<  and  can  be  expressed 


k  *  NRk  *  1  NIk 

k=0,  ...  n/2 

(4.5a) 

k  ‘  SRk  *  1  Slk 

k=0,  ...  n/2 

(4.5b) 
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Fk*M 


K  -  <sk^> 


,  k=0,  ...  n/2. 


(4.10) 


(The  -w2  term  is  included  in  the  denominator  since  z  +  s^,  *=0, 
...  n-1,  is  measured  and  its  OFT  is  -wk  Zk  +  Sk,  k=0,  ...  n/2.) 

We  wish  to  know  the  minimum  value  in  |  Fk j 2/ j Zk | 2,  but  when  noise 
is  present,  we  must  estimate  this  using  Q£.  Qk  is  a  random  vari¬ 
able  since  it  is  a  function  of  the  random  variables,  N|<  and 
Sfc.  Since  the  parameter  a0  will  be  estimated  using  the  minimum 
value  of  Q£,  the  estimator  for  a0  is  a  random  variable.  Its 
probability  distribution  is  a  function  of  the  probability  dis¬ 
tribution  of  Q£.  The  probability  distribution  of  Qk  is  quite 
complicated  because  of  the  form  of  its  functional  dependence  on 
and  S^.  However,  the  moments  of  Qk  can  be  easily  estimated. 
The  first  step  toward  accomplishing  this  task  is  to  compute  the 
moments  of  the  numerator  and  denominator  of  Q£. 

Consider  the  numerator  in  Equation  4.10.  Its  mean  and  vari¬ 
ance  can  be  found  using  the  Expressions  4.8.  Note  that  F^  has 
been  assumed  a  nonrandom  (though  unknown)  constant.  The  mean  of 
|Fk  +  NjJ2  is  the  expected  value  of  Fk  +  N^  times  its  complex 
conjugate.  This  is 

E  [|  Fk  +  Nk[2]  =  |  Fk{  2  +  nan2  At2,  k=0,  ...  n/2.  (4.11a) 

The  variance  of  | Fk  +  Nk|2  is  the  mean  square  of  ' F^  +  Nk|2  minus 
the  square  of  the  mean.  This  is 


V[|Fk  +  Nk|2]  =  n2on4  At4  +  2nan2  At2  \Fy 


+  26ko  nan2  AtZtFR0  +  FI0^’  k=0‘  **•  n/2 

(4.11b) 


where  fyo  Is  the  Kronecker  delta.  The  mean  and  variance  of  the 
denominator  of  Equation  4.10  can  also  be  found  and  these  are 

2  2 

E [j zk  -  Sk/w2)  ]  =  |Zk|  +  no2  at2/w4,  k=0,  ...  n/2  (4.12a) 

2 

V[|Zk  -  Sk/w2|  ]  =  n2as4At4/w®  +  2na2At2|Zk| 2/w4 

+  2  6k0  n<Js2At2{ZRQ  +  z j q ) /wk »  k=0,  ...  n/2 

(4.12b) 

So  far,  the  moments  established  in  Equations  4.11  and  4.12  are 
exact,  given  the  assumptions  concerning  n^  and  se,  4=0,  ... 
n-1,  listed  previously.  We  now  form  approximations  for  the 
moments  of  Qk  in  Equation  4.10. 

Recall  that  if  X  and  Y  are  independent  random  variables  and 
if  Z  is  defined  as  the  ratio  of  these,  then  the  first  two  moments 
of  Z  can  be  approximated.  Let 

Z  =  X/Y.  (a) 

Then 

E[Z]  =  E[X]/E[Y]  (b) 

V[Z]  -  ( (E[X] )2  V[Y]  +  (E[Y])2  V[X])/(E[Y])4.  (c) 

Based  on  these  formulas,  the  approximate  values  for  the  mean  and 
variance  of  Qk  can  be  written 

%  2  /  ?  2 

E[Qk]  *  E[(Fk  +  Nk!  yE[|Zk  -  Sk/w‘l  ],  k-0,1,  ...  n/2 

(4.13a) 
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and 


2  2  2  2  2  2 

!E[|Fk*M  ]}<[IV*kl  HE[|Z„-Sk/«‘|  ])  V[|VV"kl  ! 

V[QJ] - — - — - — - — — —  L- 

(E[|Zk-Sk/.k?|  ]| 

k=0,  ...  n/2. 

(4.13b) 

These  moments  characterize  the  distribution  of  Q£,  k=0,  ...  n/2, 
but  do  not  completely  define  it.  As  mentioned  above,  the  dis¬ 
tribution  of  is  quite  complicated.  For  the  purposes  of  this 
approximate  analysis,  however,  we  can  assume  that  Q£  has  a  dis¬ 
tribution  which  takes  a  relatively  simple,  common  form.  We  as¬ 
sume  Q£  to  be  normally  distributed.  Then  the  distribution  of  Q£ 
is  fully  described  by  E[Q£]  and  V[Q£].  With  this  assumption,  the 
analysis  can  proceed. 

The  estimated  value  of  the  parameter  a0  depends  on  the  mini¬ 
mum  value  of  Q£,  k=0,  ...  n/2.  Let  be  the  frequency  where 
assumes  its  minimum.  Then  we  can  write 

P(^n  *  wk)  =  P(Q£  is  the  smallest  of  the  values 
Qj  ,  j=0,  ...  n/2) 

•  p<n  <  oj) 

■  p(  n  -  Qj  <.  °>  k=°*  «/* 

j*k 

(4.14) 
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He  re  1 1 j  ^  <.  Qj  simply  refers  to  the  joint  event  where  is 

equal  to  or  less  than  each  Q‘. ,  j*k.  The  last  step  simply  rear- 

J 

ranges  the  inequality.  Equations  4.8  and  4.9  show  that  the  ran¬ 
dom  variables  Q£  and  Qj  are  independent  when  j*k.  Therefore, 
Equation  4.14  can  be  written 

=TTp«k  -oj  <«) 

j*k 

=TTF0'  -  0!(0)  k=0,  ...  n/2 


(4.15) 

where  Fgj*  _  Qj(x)  is  the  cumulative  distribution  function  of  Q|(  - 
Oj.  (Since  is  assumed  normally  distributed,  then  Q£  -  Qj  will 
also  be  normally  distributed.) 

Finally,  Equation  3.27  in  Section  3  establishes  the  formula 
for  a0  as 


ao  *  "V 


(4.16) 


In  terms  of  this  formula,  the  probability  distribution  for  the 
estimator  of  a0  is 


P(A  n  =  <*)  = 


P(£lm  =  yja/m)  ,  a  =  m k=0,  ...  n/2 


otherwise 


(4.17) 


where  A0  is  defined  as  the  random  variable  used  to  estimate  a0. 
This  formula  provides  an  approximate  probabilistic  description  of 


the  character  of  the  estimator  A0.  We  expect  that  when  the  func¬ 
tion  Q£,  k=0,  ...  n/2,  has  a  pronounced  trough  near  its  minimum, 
the  probability  distribution  of  A0  will  be  tightly  concentrated 
around  mw*.  When  Q£ ,  k=0,  ...  n/2,  is  relatively  shallow,  the 
distribution  of  A0  may  show  considerable  spread. 

In  order  to  establish  confidence  intervals  on  the  underlying 
value  of  a0,  a  probability  relation  between  the  underlying  value 
and  its  estimator,  A0,  must  be  established.  This  relation  can  be 
solved  for  a0  to  establish  confidence  intervals  on  a0.  For  the 
approximate  purposes  of  this  analysis,  we  simply  assume  that  A0 
is  a  normally-distributed  random  variable  with  mean  a0.  (The  as¬ 
sumption  that  A0  has  mean  ao  is  not  exactly  correct,  as  may  be 
noted  from  Equation  3.26.)  Then  confidence  intervals  on  the  un¬ 
derlying  constant,  a0,  can  be  written 

C(mw^  -  <Jj  <  aQ  <  +  b  =  2#(b)  -  1  b  >  0 

(4.18) 

where  o*  is  the  estimated  variance  of  A0  obtained  using  Equation 
4.17,  t[  •)  is  the  standard,  normal  cumulative  distribution  func¬ 
tion;  and  b  is  a  constant  related  to  the  confidence  level  of  the 
established  intervals.  This  concludes  the  error  analysis  on  the 
parameter  a0.  The  error  analysis  on  aj,  j=l,  ...  M+-1,  can  now 
be  performed. 

Error  analysis  on  the  variable  a!  will  be  performed  for  the 
case  where  M=0.  This  analysis  will  be  conditioned  on  our  know¬ 
ledge  of  the  value  (or  estimate)  of  ao>  This  approach  is  taken 
since  Equations  3.28,  3.29,  and  3.34  through  3.40  rely  on  this 
knowledge  of  ao.  When  noise  Is  present  in  the  structural  input 
and  response  measurements,  the  formula  used  to  estimate  a!  be¬ 
comes 
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(4.19) 


lEqk 


The  random  variable  has  replaced  Qk  in  Equation  3.29;  there¬ 

fore,  is  the  random  variable  estimator  for  the  underlying 
parameter  a^.  The  mean  and  variance  of  A*,  conditioned  on  a0, 
can  be  established  based  on  the  mean  and  variance  of  Q£  given  in 
Equations  4.13.  These  are 

e^|40]  ■  <E  wk4)  1  <E  E!qkl  "k  •  E  "(“n,2 

k  k  k 

nA?|a0].  (L«k4)'2  (E  V[Qkl»k4l 

k  k 

With  this  information,  an  approximate  conditional  probability 
distribution  for  A^  given  a0  can  be  established.  We  take  the 
conditional  distribution  of  Ag  to  be  normal  with  the  mean  and 
variance  listed  above.  Using  these  assumptions,  the  confidence 
intervals  on  a^  can  be  written 

C(E[Af|  aQ]  -  b^V[Af  |aQ]  <  a*  <  E [Af |  aQ]  ♦  |  aQ) 

=  2*(b)  -  1  b  >  0 

(4.21) 

For  various  values  of  b,  and  using  the  estimated  moments  of  Equa¬ 
tions  4.20,  Equation  4.21  sets  the  confidence  limits  for  a*. 


(4.20a) 

(4.20b) 
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We  now  perform  an  error  analysis  on  the  variables  aj  and  a2 
for  the  case  where  M=l.  As  above,  we  condition  our  analysis  on 
knowledge  of  the  value  of  a0.  Equation  3.40  establishes  an  ex¬ 
pression  for  the  two-dimensional  vector  {b}.  Its  elements,  bj 
and  b2,  are  simply  related  to  aj  and  a2>  so  we  concentrate  our 
preliminary  analysis  on  {b}.  When  the  value  of  a0  is  specified, 
the  matrix  [XL]  in  Equation  3.40  is  known.  When  Q£  replaces  Qk 
in  Equation  3.38c,  (b}  must  be  replaced  by  a  random  vector.  We 
perform  the  specified  replacement  to  obtain 

{B  }  =  [R  ]  {X^ }  (4.22) 

where  {B}  is  the  estimator  random  vector  for  {b } ,  { }  is  the 
vector  whose  form  is  Equation  3.38c  with  replaced  by  Q£,  and 

t«]  ■  ([X,]T  I*!])’1  [*i!T 

where  [Xx]  is  defined  in  Equation  3.38b.  The  conditioned  mean 
value  of  (B }  is 

E[{B}|a0]  =  [R]  E(X'}  (4.23) 

where 

'(OkjJ  -  mW  -  “k/ 

EK,»il  -  *ZK2  ■  "k2^)2 
« 

Ef»i2)  -  ”2<»„2  -  “2/ 

The  expression  for  E[Q^  ]  is  given  in  Equation  4.13a.  Equation 
4.23  provides  the  mean  value  of  Bj  and  B2,  given  a0. 
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The  conditioned  covariance  matrix  for  {B}  is  obtained  using 
the  following  procedure:  1)  Find  the  mean  value  of  the  product 
{8}  {B  }T.  This  is  done  using  Equations  4.22,  3.38c  and  4.13, 
and  by  noting  that  Q£  and  Qj  are  independent  for  k*j.  2)  Sub¬ 
tract  from  this  the  product  E[ {B } J  a0]  (E[{B}|  a 0 ] )T.  The  re¬ 
sult  is 


[KBB|a0]  ’  ("l  W  Ir)T 


(4.24) 


where  the  diagonal  elements  of  [K]  are  the  variances  of  the  ran 
dom  variables  ,  k=klt  ...  k2,  and  the  off-diagonal  terms  are 
zero. 


V[Q£  ]  0  ...  0 

0  V[Q‘i+1]  ...  0 


(4.24a) 


0  0  V[Q'  ] 

k2  J 


The  diagonal  elements  in  [ Kgg i a Q ]  are  the  variances  of  Bj  and 

B2,  given  a0.  Equations  4.23  and  4.24  define  the  first-  and 

second-order  moments  of  {B }.  These  can  be  used  to  establish  the 

moments  for  the  estimators  of  aj  and  a2. 

Let  A.  be  the  random  variable  of  the  estimator  a-,  j  =  1 , 2 . 

J  j 

Then,  according  to  Equation  3.41,  we  can  write 


A*  =  Bj,  A*  -  b|/Bj  (4.25) 

The  conditional  moments  of  A^  can  be  obtained  by  directly  taking 
the  mean  and  variance  on  both  sides  of  the  first  of  Equations 
4.25. 


/ 

h 


3 


42 


(4.26a) 


E[»l  l*o)  ■  £(B1  iao! 
v  E*f  I  a0  ]  *  V[Bj  |a„J 


(4.26b) 


The  moments  on  the  right-hand  side  are  defined  in  Equations  4.23 


and  4.24.  The  moments  of  A^  must  be  approximated  since  the  ex¬ 
pression  for  A|  involves  a  ratio  and  a  power  of  Bj  and  B2.  Usii 
Equations  a,  b  and  c,  listed  previously  in  this  section,  we  can 
show  that  the  conditional  moments  of  A^  are  approximately 

E [A|  |  ] 

a*  E[B|  |  aoJ/EtB!  I  ae] 

(4.27a) 

V[A|  la0] 

2  2 
0.  (E[b|  |  a0])  V[B1  |a0]  +  (E[B1  |a0j) 

V(B||a„) 

*3* 

o 

H 

CO 

Ul 

1 

(4.27b) 

where 

E(Bil  ao  1 

3  V[B2|  a0]  +  (E[B2|  a0])2 

(4.27c) 

*[b||  aol 

-  3(V[B2)  aQ])2 

(4.27d) 

The  moments  in  Equations  4.27c  and  4.27d  can  be  obtained  directly 
from  Equations  4.23  and  4.24.  These  can  be  used  to  establish 
Equations  4.27a  and  4.27b.  Equation  4.27d  uses  the  normal  ap¬ 
proximation  to  estimate  the  variance  of  B|. 

When  we  establish  the  relation  between  each  underlying  pa¬ 
rameter,  aj,  and  its  estimator  random  variable,  A j ,  j*l,2, 
and  when  we  find  the  probability  distribution  of  the  Aj ,  j=l,2, 
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we  can  write  confidence  intervals  on  the  aj.  For  the  approxi¬ 
mate  purposes  of  the  present  analysis,  we  simply  assume  that  the 
A?,  j=l,2,  are  normally-distributed  random  variables.  Then  the 

J 

confidence  intervals  on  the  aj ,  j=l,2,  can  be  written 

C(E[A^  |  a0]  -  C Vv[a|  I  a0]  <  aj  <  E[a|  |  aQ]  +  \  aQ]) 

=  2«(c)  -  1,  j  =  1 , 2 ,  c  >  0  (4.28) 

This  formula  reflects  our  confidence  that  a?  falls  within 
_+  c^V [A|  |  a  0]  of  its  estimated  value. 

Equations  4.18,  4.21,  and  4.28  establish  confidence  inter¬ 
vals  on  the  parameters  aj,  j=0,l,2.  The  confidence  intervals 
on  ax  and  a2  are  actually  conditioned  on  the  value  of  a0.  These 
could  be  made  unconditional  by  deriving  them  from  the  joint  prob¬ 
ability  distribution  for  A0,  Alt  and  A2.  This  joint  distribution 
could  be  obtained  by  multiplying  the  conditional  distribution  of 
A  j  and  A2  by  the  distribution  of  A0.  But  this  was  not  done  in 
the  present  analysis  for  the  following  reason.  The  separation  of 
the  error  analysis  into  two  sequential  parts  follows  the  approach 
used  in  the  parameter  identification  of  Section  3.  First,  the 
value  of  a0  is  estimated.  Then  this  value  is  used  in  the  identi¬ 
fication  procedure  to  find  a*  and  a2. 

Simplifying  and  approximating  assumptions  regarding  the  dis¬ 
tributions  of  the  parameter  estimators,  A0,  At,  and  A2,  have  been 
made  freely  throughout  this  analysis.  These  assumptions  cer¬ 
tainly  affect  the  accuracy  of  this  error  analysis,  but  we  feel 
that  useful  results  can  still  be  obtained  using  the  formulas  pre¬ 
sented  here.  In  particular,  when  the  variance  of  the  noise  on 
the  measured  structural  input  and  response  are  small,  we  expect 
that  the  analyses  of  this  section  will  be  accurate. 

The  error  analysis  for  the  time  domain  approaches  to  param¬ 
eter  identification  are  not  included  in  this  report  for  two  rea¬ 
sons.  1)  The  time  domain  analyses  developed  in  this  report  yield 


results  which  are  more  sensitive  to  the  presence  of  noise  than 
the  frequency  domain  approaches.  Therefore,  these  approaches 
will  be  less  useful  to  us  in  practice.  2)  The  time  domain  error 
analyses  are  quite  difficult  to  implement  in  a  practical  computa¬ 
tion  scheme. 

5.0  Numerical  Examples 

In  this  section,  numerical  examples  are  presented  which  dem¬ 
onstrate  the  use  of  the  analytic  procedures  developed  in  the  pre¬ 
vious  sections.  Two  examples  demonstrating  the  time  domain  ap¬ 
proach  to  parameter  identification  are  summarized.  One  example 
showing  the  frequency  domain  approach  is  presented.  The  examples 
demonstrate  the  identification  of  the  model  parameters  for  linear 
and  hysteretic,  single-degree-of-freedom  (SDF)  structures,  both 
when  measurement  noise  is  and  is  not  present. 

The  input  used  to  excite  the  SDF  systems  in  all  the  numeri¬ 
cal  examples  is  a  decaying  exponential,  oscillatory  function.  It 
is  generated  using  the  formula 

t  N 

f(t)  =  e'a  c j  cos (Wjt  -  4> j )  ] »  0  <  t  <  T.  (5.1) 

j  =  l 


a,  c • ,  j=l,  ...  N,  and  w, ,  j=l,  ...  N,  are  constants.  <* . ,  j=l, 

J  J  J 


...  N,  are  phase  angles  which  are  random  variable  realizations; 
the  underlying  random  variables  are  independent  and  uniformly 
distributed  on  the  interval  (0,  2n).  a  is  the  input  decay  rate. 
The  cj,  j=l,  ...  N,  are  constants  determining  the  amplitudes  of 
input  components.  In  all  cases  described  below,  all  the  values 
of  the  cj  are  taken  as  equal  to  c.  w j ,  j=l,  ...  N,  are  the 
frequencies  where  the  input  components  are  generated.  In  all  ex¬ 
amples  shown  here,  the  wj,  j=l,  ...  N,  are  equally  spaced  in  an 
interval  including  the  characteristic  frequency  of  the  system  be¬ 
ing  analyzed.  The  forcing  function  defined  above  was  generated 
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at  discrete  times  for  use  in  the  numerical  examples.  Specifi¬ 
cally,  f(t)  was  evaluated  at  the  times  t  =  tA  =  ut,  i=0,  ... 
n-1.  A  computer  program,  named  FORCE,  which  generates  the  inputs 
of  Equation  5.1  was  used  in  these  numerical  examples. 

It  was  necessary  to  compute  both  linear  and  hysteretic,  SDF 
structure  responses  in  the  examples.  This  was  done  using  an  in¬ 
cremental  equations  approach,  described,  for  example,  in  Refer¬ 
ence  8.  A  computer  program,  named  BILIN,  was  written  to  execute 
the  structural  response  computations.  BILIN  can  be  used  to  find 
the  displacement,  velocity,  and  acceleration  response  of  a  bi¬ 
linear  hysteretic  structure  to  an  arbitrary  input.  It  also  com¬ 
putes  the  energy  dissipated  by  the  structure  during  its  re¬ 
sponse. 

Two  basic  approaches  for  solving  the  parameter  identifica¬ 
tion  problem  were  proposed  in  this  study.  These  are  the  time  do¬ 
main  and  frequency  domain  approaches.  Two  computer  programs  were 
written  to  execute  the  parameter  identifications.  These  are 
named  TIMEID  and  FREQID.  TIME  ID  performs  time  domain  parameter 
identification.  It  accepts  an  input  signal  from  FORCE  and  a  com¬ 
puted  response  from  BILIN  and  then  does  one  of  two  things.  It 
either  submits  the  data  for  an  immediate  parameter  identifica¬ 
tion,  or  it  adds  white  noise  to  the  input  and/or  output,  and  then 
submits  the  data  for  parameter  identification.  The  program  can 
perform  the  parameter  identification  directly  using  Equations  3.4 
through  3.12,  or  it  can  perform  the  identification  using  the  in¬ 
tegrated  equations  approach,  Equations  3.13  through  3.18.  The 
program  can  identify  the  parameters  of  the  second-order  model 
(M=0)  or  the  third-order  model  (M=l)  of  Equation  3.4. 

The  computer  program  FREQID  performs  frequency  domain  param¬ 
eter  identification.  It  accepts  an  input  signal  from  FORCE  and  ? 
computed  response  from  BILIN.  If  desired,  FREQID  next  adds  white 
noise  signals  to  the  input  and/or  response.  Then  it  performs  the 
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necessary  Fourier  transforms  and  other  data  operations.  Follow¬ 
ing  this,  the  parameter  identification  is  executed.  Equations 
3.19  through  3.41  are  used.  The  program  can  identify  the  param¬ 
eters  of  the  second-order  model  (M=0)  or  the  third-order  model 
(M=l)  of  Equation  3.4. 

Once  the  parameters  of  Equation  3.4  have  been  estimated,  the 
energy  dissipated  by  the  model  is  computed;  and  this  is  compared 
to  the  energy  dissipated  by  the  actual  system,  as  computed  in 
B IL IN.  This  computation  is  the  one  discussed  in  Section  2  and 
given  by  Equation  2.5b.  The  energy  computations  are  performed  in 
programs  ENED10,  for  second-order  systems,  and  ENED13,  for  third- 
order  systems.  The  computations  are  performed  using  an  incremen¬ 
tal  form  of  the  governing  Equation  3.4. 

The  notations  for  the  parameters  used  in  specifying  the  nu¬ 
merical  examples  are  those  used  in  the  text.  Some  additional 
notations  are  defined  here,  c  is  the  viscous  damping  in  an  SDF 
system,  k  is  the  initial  stiffness  in  a  bilinear  hysteretic 
structure.  ky  is  the  yield  stiffness  of  a  bilinear  hysteretic 
structure.  zy  is  the  yield  displacement  of  a  bilinear  hyster¬ 
etic  structure.  zmax  is  the  maximum  displacement  of  an  SDF 
structure. 

In  the  following  numerical  examples,  four  basic  problem 
types  are  solved.  These  are  summarized  below. 

(1)  An  input  is  generated  using  Equation  5.1.  The  input  is  used 
to  excite  a  linear  SDF  system  with  viscous  damping.  The 
structural  input  and  response  are  stored,  and  no  noise  sig¬ 
nals  are  added  to  the  input  and  response.  Then  the  input 
and  response  are  used  to  identify  the  model  parameters, 

aj,  j=0,  ...  M+l,  from  Equation  3.4. 

(2)  An  input  is  generated  as  in  1,  above,  but  here  the  response 
of  a  bilinear  hysteretic  system  is  computed.  No  noise  sig¬ 
nals  are  added  to  the  input  and  response.  The  input  and  re¬ 
sponse  are  used  to  identify  the  model  parameters,  aj,  j=0, 
...  M+l. 
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(3)  This  case  is  the  same  as  1,  above,  but  here  noise  signals 
are  added  to  the  input  and  response.  Then  the  parameters 
aj,  j=0,  ...  M+l  are  identified. 

(4)  This  case  is  the  same  as  2,  above,  but  here  noise  signals 
are  added  to  the  input  and  response.  Then  the  parameters 
aj,  j=0,  ...  M+l  are  identified. 

In  all  four  cases,  the  parameters  of  the  second-order  model,  a0 
and  (M=0),  and  the  third-order  model,  a0,  alt  and  a2  (M=l), 
are  identified.  Moreover,  the  energy  dissipated  by  the  model 
system  with  the  identified  parameters  is  computed. 

5.1  Example  1 

This  example  solves  a  parameter  identification  problem  using 
the  direct  time  domain  approach,  summarized  in  Equations  3.4 
through  3.12.  The  parameters  of  the  shock  input  are  listed  in 
Table  5.1.  A  typical  forcing  function  history  generated  using 
these  parameters  is  shown  in  Figure  5.1a.  The  derivative  of  the 
forcing  function  is  shown  in  Figure  5.1b. 

The  responses  of  some  SOF  systems  to  the  shock  input  in  Fig¬ 
ure  5.1a  were  computed.  First,  the  response  of  a  linear  system 
was  computed  for  analysis  in  cases  1  and  5.  The  computed  dis¬ 
placement  response  is  plotted  versus  time  in  Figure  5.2a.  The 
SDF  structure  spring  restoring  force  versus  displacement  is 
plotted  in  Figure  5.2b.  The  very  slightly  nonlinear  response  of 
an  SDF  structure  was  computed  for  analysis  in  cases  2  and  6,  but 
this  response  is  not  shown.  A  more  nonlinear  response  was  com¬ 
puted  for  analysis  in  cases  3  and  7.  The  displacement  response 
versus  time  is  plotted  in  Figure  5.3a,  and  the  spring  restoring 
force  versus  displacement  is  shown  in  Figure  5.3b.  The  first 
figure  shows  a  residual  plastic  displacement  as  the  response  vi¬ 
brations  diminish.  The  second  graph  shows  the  permanent  set  as 
lateral  displacement  of  the  horizontal  axis  intercept.  Finally, 
a  very  nonlinear  response  was  computed  for  analysis  in  cases  4 
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FORCE/SEC  o  .  FORCE  -  KIPS 


4000.0 


Figure  5.1a.  Decaying  exponential  forcing  function 
for  Example  1 


0.0  0.5  1.0  1.5 

TIME  -  SEC 

Figure  5.1b.  Derivative  of  forcing  function  in 
Figure  5.1a 
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DISPLACEMENT 


Table  5.1.  Parameters  and  Results  for  Example  1 


Input  (Structural  Excitation) 


Case  Number. 

|U  _  W,  1 

N  a  C  1 

*N 

n  At 

1  through  8 

50  1.0  300  32.4  1256.4 

256  0.005 

Structure  Parameters 

Case  Number 

m  c  k  ^y 

Zy 

E0 

1,5 

0.259  7.77  5829  - 

00 

23292 

2,6 

0.259  7.77  5829  0 

7700 

23588 

3,7 

0.259  7.77  5829  '  0 

7000 

23674 

4,8 

0.259  7.77  5829  0 

5000 

21959 

Identification  Parameters  and  Results 

Case  Number 

M  a0  al  a2 

E0 

1 

2  5829  7.77 

21533 

2 

2  5825  7.79 

21533 

3 

2  4574  8.13 

21779 

4 

2  1620  9.69 

23844 

5 

3  5829  7.77  0.0 

21533 

6 

3  5823  7.85  0.0 

21498 

7 

3  2527  99.65  0.0158 

20824 

8 

3  610  108.11  0.0175 

18666 

and  8.  The  displacement  response  versus  time  is  shown  in  Figure 
5.4a  and  th*  spring  restoring  force  versus  displacement  is  shown 
in  Figure  5.4b.  A  considerable  permanent  set  is  evident,  as  the 
motion  diminishes,  from  the  first  figure.  The  second  figure 
shows  that  large  plastic  deformations  occur  in  the  structure  in 
both  directions  of  motion. 

The  energy  dissipated  by  each  structural  system  is  listed 
with  the  structural  parameters  in  Table  5.1.  Eg  is  that  energy 
dissipated  due  to  the  action  of  the  inelastic  spring  and  the  ac¬ 
tion  of  the  viscous  damper. 
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Using  the  forcing  function  Input,  described  above,  and  the 
computed  responses,  the  parameters  of  the  structural  systems  were 
identified.  The  results  of  the  parameter  Identifications  are 
given  in  Table  5.1.  The  energy  dissipated  when  the  identified 
systems  respond  to  the  shock  input  is  listed  in  Table  5.1  next  to 
the  identified  parameters. 

Figures  5.5a  through  5.9b  show  the  computed  responses  of 
some  of  the  Identified  systems.  The  figure  titles  indicate  which 
systems  generate  the  responses  shown.  The  top  (or  "a")  figures 
show  the  computed  responses  versus  time.  The  bottom  (or  "b") 
figures  show  the  computed  restoring  forces,  spring  force  plus 
damper  force,  versus  displacements. 

Figure  5.10  compares  three  responses.  These  are:  1)  the 
actual  nonlinear  structural  response  obtained  using  BILIN  in  the 
slightly  nonlinear  cases  2  and  6,  2)  the  response  executed  by  the 
identified  model  described  in  case  2,  and  3)  the  response  exe¬ 
cuted  by  the  identified  model  described  in  case  6.  The  two  lat¬ 
ter  responses  practically  overlay  and  very  nearly  equal  the  ac¬ 
tual  response. 

Figure  5.11a  compares  the  actual  nonlinear  response  of  cases 
3  and  7  to  the  second-order  (M=0)  model  response  of  case  3.  Fig¬ 
ure  5.11b  compares  the  actual  nonlinear  response  of  cases  3  and  7 
to  the  third-order  (M=l)  model  response  of  case  7.  Both  models 
simulate  response  amplitudes  quite  well;  the  third-order  model  is 
slightly  closer  than  the  second-order  model  in  that  the  third- 
order  model  provides  a  slightly  better  phase  match  to  the  actual 
response  than  the  second-order  model. 

In  this  numerical  example,  no  cases  are  included  where  noise 
was  added  to  the  forcing  function  input  and/or  the  acceleration 
response.  Such  cases  were  analyzed,  but  the  results  were  so  poor 
that  they  are  not  summarized  here.  These  results  showed  that  the 
direct,  time  domain  parameter  identification  technique  is  not 
effective  in  parameter  analysis  when  recording  noise  is  present. 
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Figure  5.10.  The  comparison  between  measured  response  (light 
line)  and  identified  response  (dark  line)  for 
both  second-order  (case  3)  and  third-order  (case 
approximate  systems.  (The  identifier  responses 
are  practically  the  same.) 


The  numerical  examples  summarized  here  show  that  the  direct, 
time  domain  parameter  identification  technique  can  be  used  effec¬ 
tively  when  the  measured  signals  are  noise  free.  This  is  best 
confirmed  by  reference  to  Figures  5.10,  5.11a,  and  5.11b.  These 
show  that  the  linear  model  response  can  be  made  to  match  the  non¬ 
linear  response  well. 

The  energy  dissipation  results,  summarized  in  Table  5.1, 
show  that  the  third-order  models  provide  the  best  simulation  for 
a  nonlinear  hysteretic  system. 

5.2  Example  2 

The  parameter  identification  problem  is  solved  in  this  exam¬ 
ple  using  the  time  domain  approach  with  the  integrated  equation 
of  motion.  This  approach  was  developed  in  Equations  3.13  through 
3.18.  The  forcing  function  input  used  to  excite  the  systems  in 
this  example  is  the  same  one  used  in  Example  1.  The  forcing 
function  parameters  are  listed  in  Table  5.2.  In  the  present 
case,  though,  the  time  increment  where  the  forcing  function  was 
evaluated  is  At  *  5  x  10“5  sec.  It  was  found,  in  working  several 
numerical  examples,  that  a  time  increment  of  extremely  small  mag¬ 
nitude,  relative  to  the  natural  period  of  the  system  under  inves¬ 
tigation,  is  necessary  to  make  the  system  identification  work 
when  noise  is  present  in  the  inputs  and  outputs. 

The  structural  systems  analyzed  in  the  6  cases  of  the  pres¬ 
ent  numerical  example  are  summarized  in  Table  5.2.  The  response 
of  a  linear  system  was  computed  for  analysis  in  cases  1  and  4. 

The  responses  of  moderately  and  strongly  nonlinear  systems  were 
computed  for  analysis  in  cases  2  and  5  and  cases  3  and  6,  respec¬ 
tively. 

In  all  cases  in  this  numerical  example,  noise  was  added  to 
the  generated  forcing  function  input  and  the  measured  accelera¬ 
tion  response.  The  standard  deviations  of  the  band-limited, 
white  noise,  random  processes  added  to  the  inputs  and  responses 
are  given  in  Table  5.2.  Figures  5.12a  and  5.12b  show  segments  of 


Figure  5.12a.  Noisy  and  noise-free  forcing  function  signals. 
The  noisy  signal  is  used  in  Example  2. 


Figure  5.12b.  Noisy  and  noise-free  acceleration  response 

signals.  The  noisy  signal  is  used  in  Example  2 


Table  5.2.  Parameters  and  Results  for  Example  2 

Input  (Structural  Excitation) 

Case  Number  N  a  c  W1  WN  n  At 

1  through  6  50  1.0  300  32.4  1256.4  256  5xl0~5 

Structural  Parameters 


Case  Number 

m 

_c _ 

k 

Z ' 

_^D_ 

1,4 

0.259 

7.77 

5829 

09 

22.7 

2,5 

0.259 

7.77 

5829 

0 

2000 

22.7 

3,6 

0.259 

7.77 

5829 

0 

1000 

22.7 

Noise  Parameters 

Case  Number  gf  CTs 

1  through  6  150  475 

Identification  Parameters  and  Results 


Case  Number 

M 

Jo_ 

al 

a2 

ed 

1 

2 

5809 

7.89 

22.8 

2 

2 

4526 

14.20 

26.5 

3 

2 

1563 

18.19 

24.6 

4 

3 

5737 

9.11 

0 

24.1 

5 

3 

3851 

70.79 

0.00141 

28.9 

6 

3 

1056 

18.50 

0.00147 

the  generated  forcing  function  and  the  computed  acceleration 
response  with  noise  added.  A  typical  segment  of  the  white  noise 
realization  added  to  the  input  and  response  is  shown,  alone,  in 
Figure  5.12c. 

Parameter  identifications  were  executed  using  the  noisy  in¬ 
puts  and  responses  summarized  above.  The  results  are  listed  at 
the  bottom  of  Table  5.2.  Comparison  of  the  actual  values  for 
energy  dissipated  and  the  identified  values  for  energy  dissi¬ 
pated,  in  Table  5.2,  shows  that  the  identified  models  provide  a 
reasonable  estimate  for  energy  dissipated  in  the  actual  system. 
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The  capability  to  establish  such  an  estimate  in  the  presence  of 
noise  provides  an  Improvement  over  the  parameter  identification 
method  used  in  Example  1. 

5.3  Example  3 

In  this  example,  a  parameter  identification  problem  is 
solved  using  the  frequency  domain  approach,  developed  in  Equa¬ 
tions  3.19  through  3.41.  The  shock  input  used  in  this  example  is 
a  decaying  exponential,  oscillatory  input.  The  input  parameters 
are  listed  in  Table  5.3.  These  parameters  are  the  same  for  all 
eight  cases  worked  in  this  example.  A  forcing  function  generated 
using  the  Table  5.3  parameters  is  shown  in  Figure  5.13a.  Param¬ 
eter  identifications  were  performed  both  when  noise  was  and  was 
not  present  on  the  forcing  function  input.  The  standard  devia¬ 
tions  of  the  noise  random  processes  for  all  eight  cases  worked  in 
this  example  are  listed  in  Table  5.3.  Note  that  cases  2,  4,  6 
and  8  are  noisy;  cases  1,  3,  5  and  7  are  not  noisy.  The  forcing 
function  input  with  a  white  noise  superimposed  on  it  is  shown  in 
Figure  5.13b. 

Two  structural  responses  to  the  input  of  Figure  5.13a  were 
computed.  The  first  response  was  for  a  linear,  SDF  system;  the 
structure  parameters  are  listed  in  Table  5.3.  This  response  was 
used  in  example  cases  1,  2,  5  and  6.  The  second  response  was  for 
a  bilinear  hysteretic  system;  the  structural  parameters  are 
listed  In  Table  5.3.  This  response  was  used  in  example  cases  3, 
4,  7  and  8.  The  acceleration  response  of  the  nonlinear  system  is 
shown  in  Figure  5.14a.  Parameter  identifications  were  performed 
both  when  noise  was  and  was  not  present  on  the  acceleration  re¬ 
sponse.  The  standard  deviations  of  the  white  noise  random  pro¬ 
cesses  are  listed  in  Table  5.3  for  all  8  cases.  Figure  5.14b 
shows  the  nonlinear  acceleration  response  from  Figure  5.14a  with 
a  measurement  white  noise  superimposed.  Figure  5.15  shows  the 
graph  of  restoring  force  (spring  plus  damper  forces)  versus  dis¬ 


placement. 


Table  5.3.  Parameters  and  Results  for  Example  3 


Input  (Structural  Excitation) 


Case  Number 

N_ 

a 

c 

n 

At 

1  through  8 

50 

0.025 

1.0 

0.8  1.2 

1024 

0.3 

Structure  Parameters 

Case  Number 

m 

c 

k  ky  zy 

zmax 

E0 

1,2 

1.0 

0.2 

1.0 

—  00 

12.5 

513.1 

3,4 

1.0 

0.2 

1.0  0 

.5  10.0 

12.7 

532.6 

5,6 

1.0 

0.2 

1.0 

—  00 

12.5 

513.1 

7,8 

1.0 

0.2 

1.0  0 

.5  10.0 

12.7 

532.6 

Noise  Parameters 

Case  Number 

aJ_ 

^s_ 

1 

0.0 

0.0 

2 

0.2 

0.5 

3 

0.0 

0.0 

4 

0.2 

0.5 

5 

0.0 

0.0 

6 

0.2 

0.5 

7 

0.0 

0.0 

8 

0.2 

0.5 

Identification  Parameters  and 

Results 

Case  Number 

M 

a0 

al 

a2 

E0 

1 

2 

.9458 

.1902 

544.8 

2 

2 

1.010 

.2660 

452.5 

3 

2 

.9214 

.1617 

590.2 

4 

2 

.9074 

.1374 

644.1 

5 

3 

.9458 

.3116 

.08083 

501.9 

6 

3 

.9270 

.2825 

.06461 

519.6 

7 

3 

.9214 

.2293 

.02744 

543.1 

8 

3 

.9074 

.2196 

.02430 

554.6 
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Figure  5.13a.  Signal  used  to  simulate  the  actual  forcing 
function  for  an  SDF  system  in  Example  3 


Figure  5.13b.  Signal  used  to  simulate  the  measured  forcing 
function  for  an  SDF  system.  (Includes  noise.) 
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300 


Figure  5.14b.  Signal  used  to  simulate  measured  acceleration 
response.  Signal  of  Figure  14a  plus  noise. 
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RESTORING  FORCE  -  KIPS 


Figure  5.15.  Total  restoring  force  of  nonlinear  system  versus 
displacement.  (Cases  3,  4,  7  and  8) 
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The  energy  dissipated  by  the  structure  in  each  case  of  the 
numerical  example  is  listed  in  Table  5.3  under  the  heading 
"Structural  Parameters." 

Parameter  identifications  were  performed  for  the  eight  cases 
described  above.  The  identified  parameters  are  listed  at  the 
bottom  of  Table  5.3.  The  energy  dissipated  by  each  of  the  iden¬ 
tification  systems  is  listed  next  to  the  identified  parameters. 

Figures  5.16a  and  5.16b  show  comparisons  of  the  responses 
computed  from  the  identified  linear  models  to  the  nonlinear 
response  of  cases  3,  4,  7  and  8.  Figure  5.16a  compares  the 
second-order  model  (M=0)  response  to  the  actual  response.  Figure 
5.16a  compares  the  third-order  model  (M=l)  response  to  the  actual 
response.  The  plots  show  that  the  third-order  model  provides  a 
better  visual  match  to  the  actual  response  than  the  second-order 
model.  However,  both  linear  models  provide  a  good  simulation  of 
the  nonlinear  response. 

Figures  5.16a  and  5.16b  make  it  apparent  that  the  frequency 
domain  parameter  identification  technique  provides  acceptable 
results,  at  least  in  those  cases  considered  here. 

The  results  show  that  the  third-order  model  is  preferable  to 
the  second-order  model  for  estimation  of  the  energy  dissipated  by 
the  nonlinear  response. 

6.0  Discussion  and  Conclusions 

The  main  objective  of  this  investigation  was  to  establish  a 
method  for  measuring  the  damage  accunulation  in  structures.  We 
assumed  that  structural  damage  is  related  to  the  energy  dissi¬ 
pated  during  a  strong  motion  event  and  attempted  to  identify  this 
quantity. 

The  estimate  for  energy  dissipated  by  a  structure  was  ob¬ 
tained  by  identifying  the  parameters  for  a  model  of  the  damaged 
structure,  then  using  the  model  to  compute  the  energy  dissipated 
during  the  strong  motion  event.  The  model  chosen  to  simulate  the 
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response  of  the  damaged  structure  is  a  system  of  linear  differen¬ 
tial  equations.  Specifically,  second-  and  third-order,  linear, 
ordinary,  differential  equations  were  used  to  model  the  motions 
of  single-degree-of-freedom  (SDF)  hysteretic  systems. 

The  parameters  of  the  model  were  estimated  using  three  ap¬ 
proaches.  Two  of  these  approaches  are  time  domain  least  squares 
approaches  where  the  modeling  error  is  minimized  with  respect  to 
the  measured  data.  The  other  is  a  frequency  domain  least  squares 
approach. 

Many  numerical  examples  were  solved,  and  some  of  them  are 
summarized  in  this  report.  Experience  obtained  in  solving  the 
numerical  examples  leads  to  the  following  conclusions. 

(1)  Linear  and  nonlinear  hysteretic  single-degree-of- 
freedom  (SDF)  system  responses  can  be  accurately  mod¬ 
eled  using  second-  and  third-order  linear  equations. 

(2)  A  direct,  time  domain  approach  can  be  used  to  identify 
model  parameters  when  the  force  input  and  acceleration 
response  measurements  are  not  noisy. 

(3)  A  time  domain  approach  using  integrated  equations  of 
motion  can  be  used  to  identify  model  parameters  when 
the  input  and  response  measurements  are  noisy.  How¬ 
ever,  an  accuracy  requirement  forcing  At  to  be  very 
small  when  this  approach  is  used  reduces  the  effective¬ 
ness  of  this  method. 

(4)  A  frequency  domain  approach  to  the  identification  of 
SDF  system  parameters  is  most  effective  in  the  present 
application.  It  C2n  be  used  when  measurement  noise  is 
present  in  the  input  and  response  and  leads  to  accurate 
results. 

(5)  The  third-order,  linear  model  is  preferable  to  the 
second-order  model  in  at  least  two  respects.  1)  The 
third-order  model  provides  a  better  match  to  the  non¬ 
linear  SDF  system  response  than  the  second-order 
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model.  2)  The  third-order  model  provides  a  better  es¬ 
timate  of  energy  dissipated  by  the  nonlinear  system 
than  the  second-order  model. 

The  present  approach  to  frequency  domain  parameter  identifi 
cation  and  estimation  of  energy  dissipated  by  a  hysteretic  SDF 
system  yields  accurate  results.  However,  the  degree  of  accuracy 
in  the  present  analysis  can  be  improved.  In  particular,  the 
present  analytical  approach  would  be  improved  by  introduction  of 
a  numerical  search  procedure  for  location  of  the  system  param¬ 
eters  which  correspond  to  a  least  squares  solution.  Such  a 
change  in  the  analysis  would  eliminate  the  need  for  the  simplify 
ing  assumptions  used  in  this  study. 

Future  investigations  should  attempt  to  improve  the  models 
established  in  this  study  and  consider  more  complicated  systems. 
The  models  considered  in  this  study  can  be  made  more  elaborate 
and,  possibly,  more  accurate.  For  example,  linear  models  with 
time  varying  parameters  can  be  considered;  and  nonlinear  models 
can  be  considered.  (The  drawback  in  the  use  of  nonlinear  models 
is  that  the  frequency  domain  analysis  can  no  longer  be  used.) 

The  analysis  of  energy  dissipated  in  two-degree-of-f reedom  struc 
tures  should  be  pursued. 

Future  studies  should  include  numerical  examples  of  error 
analysis.  The  present  investigation  derived  the  equations  for 
error  analysis,  but  no  numerical  computations  of  confidence  in¬ 
tervals  were  included  in  the  examples. 
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